C     =============================================================70       
C     This subroutine is to solve the heat and mass transfer in DCMD,
C     based on the solution of differential macroscopic balance for 
C     nonisothermal system (referred to B.R.Bird, Transport 
C     Phenomina, 2nd Ed., Page 460)
      SUBROUTINE DCMD_SOLV(OPT, LENLOC, WINLOC, TINLOC, PINLOC,
     &                    VINLOC, WOUTLOC, TOUTLOC, POUTLOC, VOUTLOC,
     &                    TWLOC, JHLOC, JMLOC, TPCLOC, EFFLOC)

C     Declare the global variables
      USE DCMD_CMOD

C     Input parameters
C     OPT - Options for flowing direction: 1 - concurrent, 
C                                          2 - countercurrent
C     LENLOC - Length of each part
C     WINLOC - Inlet mass flowrates of tube side and shell side 
C              [kg/s]
C     TINLOC - Inlet temperatures of tube side and shell side [K]
C     PINLOC - Inlet pressures of tube side and shell side [Pa]
C     VINLOC - Inlet velocity of tube side and shell side [m/s]
      INTEGER OPT
      REAL*8 LENLOC, WINLOC(2), TINLOC(2), PINLOC(2), VINLOC(2)
C
C     Outputs
C     WOUTLOC - Outlet mass flowrates of tube side and shell side
C               [kg/s]
C     TOUTLOC - Outlet temperature of tube side and shell side [K]
C     POUTLOC - Outlet pressures of tube side and shell side [Pa]
C     VOUTLOC - Outlet velocity of tube side and shell side [m/s]
C     TWLOC - Wall temperatures at the calculated node [K]
C     JHLOC - Heat transfer flux at the calculated node [J/m2]
C     JMLOC - Mass transfer flux at the calculated node [kg/m2]
C     TPCLOC - Temperature polarization coefficient
C     EFFLOC - Thermal efficient at the calculated node
      REAL*8 WOUTLOC(2), TOUTLOC(2), POUTLOC(2), VOUTLOC(2),
     &       TWLOC(2), JHLOC, JMLOC, TPCLOC, EFFLOC
C
C     Intermediate variables
      INTEGER ISIDE, I
      REAL*8 TM, dP(2), ALOC, CP(2), CSA(2), RHO(2)
      REAL*8 C, KM, THICK
      REAL*8 RE(2), NU(2), TW(2)
C      
      REAL*8 PI, HALF, ONE, TWO
      DATA PI /3.1415926D0/, HALF /5.D-1/, ONE /1.D0/, TWO /2.D0/
C
C     Functions
      REAL*8 GradSP, GradDP, DHVAP, SVP
C
C     Declare the parameters needed in ODE solver
      EXTERNAL F
      INTEGER*4 NEQN, IFLAG
      PARAMETER (NEQN = 2)
      REAL*8 T(NEQN), Z, ZEND, RELERR, ABSERR, WORK(100+21*NEQN)
      INTEGER*4 IWORK(5)     

C     ODE solver initialization
      RELERR = 1.D-4
      ABSERR = 1.D-7
      IFLAG = 1
      Z = 0.0
      ZEND = LENLOC

C     Get the initial value of integration
      DO I = 1, 2
        T(I) = TINLOC(I)
      END DO
C     Invoke ODE solver
      CALL ODE(F, NEQN, T, Z, ZEND, RELERR, ABSERR,
     &         IFLAG, WORK, IWORK)
      IF (IFLAG .EQ. 2) THEN
        DO I = 1, 2
          TOUTLOC(I) = T(I)
        END DO
      END IF

C     Retrieve data from common data block
      C = COM_PERM
      KM = COM_COND
      THICK = COM_THK
      DO ISIDE = 1, 2
        CP(ISIDE) = COM_CP(ISIDE)
        RHO(ISIDE) = COM_RHO(ISIDE)
      END DO
      CSA(1) = PI*(COM_ID(1)/TWO)**TWO
      CSA(2) = COM_CSA(2)
      dP(1) = GradDP(COM_DEQ(1), COM_RHO(1), VINLOC(1), COM_MU(1))
      dP(2) = GradDP(COM_DEQ(2), COM_RHO(2), VINLOC(2), COM_MU(2))

C     Following variables in common block will be used to 
C     data exchange with ODE functions in its first invoking
      DO ISIDE = 1, 2
        COM_WLOC(ISIDE) = WINLOC(ISIDE)
        COM_TLOC(ISIDE) = TINLOC(ISIDE)
        COM_PLOC(ISIDE) = PINLOC(ISIDE)
        COM_VLOC(ISIDE) = VINLOC(ISIDE)
      END DO

C     Compute the heat and mass transfer fluxes
C     based on inlet properties
      CALL DCMD_FLUX(WINLOC, TINLOC, PINLOC, VINLOC, 
     &              RE, NU, TW, JHLOC, JMLOC, TPCLOC, EFFLOC,
     &              IFLAG)
C     Compute the change of mass flowrate      
      ALOC = LENLOC*PI*COM_ID(1)*COM_NUM
C     Tube side
      ISIDE = 1
      TWLOC(ISIDE) = TW(ISIDE)
      WOUTLOC(ISIDE) = WINLOC(ISIDE)-JMLOC*ALOC
      VOUTLOC(ISIDE) = WINLOC(ISIDE)*VINLOC(ISIDE)/WOUTLOC(ISIDE)
      POUTLOC(ISIDE) = PINLOC(ISIDE)+LENLOC*dP(ISIDE)
      
C     Shell side
      ISIDE = 2
      TWLOC(ISIDE) = TW(ISIDE)
      WOUTLOC(ISIDE) = WINLOC(ISIDE)+JMLOC*ALOC
      VOUTLOC(ISIDE) = WINLOC(ISIDE)*VINLOC(ISIDE)/WOUTLOC(ISIDE)
      POUTLOC(ISIDE) = PINLOC(ISIDE)+LENLOC*dP(ISIDE)
      
C     Output results

      RETURN
            
      END SUBROUTINE     

C     ==================================================================       
C     Right hand side of heat transfer ODE equation in concurrent
      SUBROUTINE F(z, T, dT)
      
      USE DCMD_CMOD
      
      IMPLICIT NONE
      
      DOUBLE PRECISION z, T(*), dT(*)
! fluid properties: rho - density [kg/m3], 
!                    cp - heat capacity [J/kg-K], 
!                     k - thermal conductivity [W/m-K],
!                    mu - viscosity [Pa-s]

      REAL*8 rho(2), cp(2), k(2), mu(2)
! fluid inlet velocity [m/s]
      REAL*8 vin(2)
! membrane characteristics: d(1) - HF tube inner diameter [m], 
!                           d(2) - HF tube outer diameter [m],
!                          de(1) - tube-side equivalent diameter [m],
!                          de(2) - shell-side equivalent diameter [m],
!                         csa(1) - tube-side cross section area [m2],
!                         csa(2) - shell-side cross section area [m2],
!                              l - length of membrane [m],
!                              c - permeability [kg/m2-s-Pa],
!                             km - thermal conductivity [W/m2-K],
!                            num - number of hollow fiber tubes,
!                          thick - membrane thickness [m]
      REAL*8 d(2), de(2), csa(2), l, km, c, num, thick
! temperature polarization coefficient
      REAL*8 tau
! function declaration
!     GradDP - correlate gradient of pressure drop [Pa] * 
!     GradVP - correlate gradient of water vapor pressure [Pa/K] *
!     DHVAP - latent heat [J/kg]
      REAL*8 GradDP, GradSP, DHVAP
!     DRE - calculate Reynolds number
!     DPR - calculate Prandlt number
!     DGR - calculate Graetz number
!     DNU - calculate Nusselt number
      REAL*8 DRE, DPR, DGR, DNU
! local variables            
      REAL*8 G(2), Re(2), Pr(2), Gr(2), VisCorr(2), Nu(2), HFILM(2)
      REAL*8 TM, JM, H, tmpH
      REAL*8 ACOEFF(2), BCOEFF(2)
      INTEGER i
      
      DO i = 1, 2
        rho(i) = COM_RHO(i)
        mu(i) = COM_MU(i)
        k(i) = COM_LAMBDA(i)
        cp(i) = COM_CP(i)
        d(i) = COM_ID(i)
        de(i) = COM_DEQ(i)
        csa(i) = COM_CSA(i)
        vin(i) = COM_WLOC(i)/rho(i)/csa(i)
      END DO
      l = COM_LEN
      c = COM_PERM
      km = COM_COND
      num = COM_NUM
      thick = 0.5D0*(COM_OD(1)-COM_ID(1))
                  
      DO i = 1, 2
C       Calculate flowing mass flux
        G(i) = rho(i)*vin(i)
C       Calculate Reynolds number
        Re(i) = DRE(G(i), de(i), mu(i))
C       Calculate Prandtl number
        Pr(i) = DPR(cp(i), mu(i), k(i))
C       Calculate Graetz number
        Gr(i) = DGR(Re(i), Pr(i), de(i), l)
C       Calculate viscosity correction factor near the wall
        VisCorr(i) = 1.0
C       Calculate Nusselt number
        Nu(i) = DNU(Gr(i), Re(i), Pr(i), VisCorr(i), i)       
C       Calculate film heat transfer coefficient
        HFILM(i) = Nu(i)*k(i)/de(i)
      END DO
C     Calculate averaged membrane temperature
      TM = 0.5*(T(1)+T(2))
C     Calculate overall heat transfer coefficient
      H = c*GradSP(TM)*DHVAP(TM)+km/thick
C     Calculate temperature polarization coefficient
      tmpH = 1.0/(1.0/HFILM(1)+1.0/HFILM(2))
      tau = 1.0/(1.0+H/tmpH)

      DO i = 1, 2
        ACOEFF(i) = -1.d0/(rho(i)*cp(i))*
     &         GradDP(de(i), rho(i), vin(i), mu(i))
        BCOEFF(i) = 4.d0/d(i)*H*tau/(rho(i)*cp(i)*vin(i))
      END DO
      BCOEFF(2) = BCOEFF(2)*num*d(1)/d(2)
      
      dT(1) = ACOEFF(1)-BCOEFF(1)*(T(1)-T(2))
      dT(2) = ACOEFF(2)+BCOEFF(2)*(T(1)-T(2))
      
      RETURN

      END SUBROUTINE

      